#!/usr/bin/perl
# 2013/08/15 Compare associations between two contact map

use strict;
use warnings;
use IO::File;
use Getopt::Std;
use Carp qw(croak);
$| = 0;

if(@ARGV != 4 or $ARGV[0] eq '--help'){
	die "Usage : $0 -a [association file1] -b [association file2]\n";
}

my %opt;
getopts("a:b:", \%opt);
my $FILE_aso1 = $opt{a};
my $FILE_aso2 = $opt{b};


#---------------------------------------
# correct information of aso2
#---------------------------------------
my $total_aso2 = 0;
my %exp2;
my $fh_exp2 = IO::File->new($FILE_aso2) or die "cannot open $FILE_aso2: $!";
while($_ = $fh_exp2->getline()){
	s/\r?\n//;
	my ($peak1, $loc1, $peak2, $loc2, $interaction, $pvalue, $qvalue) = split /\t/;

#	if($qvalue eq 'NA' or $qvalue >= 0.05){
#		next;
#	}


	my ($chr1, $start1, $end1) = split /,/, $loc1;
	my ($chr2, $start2, $end2) = split /,/, $loc2;

	my $center1 = $start1 + 1500;
	my $center2 = $start2 + 1500;

	push @{$exp2{$chr1}{$chr2}}, [$center1, $center2];
	$total_aso2++;
}
$fh_exp2->close();

#---------------------------------------
# comapre to aso1
#---------------------------------------
my $total_aso1 = 0;
my $common_among_aso1 = 0;
my %already_defined;	# to check duplicated
my $fh_exp1 = IO::File->new($FILE_aso1) or die "cannot open $FILE_aso1: $!";
while($_ = $fh_exp1->getline()){
	s/\r?\n//;
	my ($peak1, $loc1, $peak2, $loc2, $interaction, $pvalue, $qvalue) = split /\t/;

#	if($qvalue eq 'NA' or $qvalue >= 0.05){
#		next;
#	}


	my ($chr1, $start1, $end1) = split /,/, $loc1;
	my ($chr2, $start2, $end2) = split /,/, $loc2;

	my $center1 = $start1 + 1500;
	my $center2 = $start2 + 1500;

	if(&check_peak($chr1, $chr2, $center1, $center2) > 0){
		$common_among_aso1++;
	}
	$total_aso1++;
}
$fh_exp1->close();

my $common_among_aso2 = 0;
foreach my $key(keys %already_defined){
	$common_among_aso2++;
}


#---------------------------------------
# report result
#---------------------------------------
print "$total_aso1\t$common_among_aso1\t$total_aso2\t$common_among_aso2\n";



sub check_peak{
	my ($chr1, $chr2, $source1, $source2) = @_;
	my @array = @{$exp2{$chr1}{$chr2}};
	my @newArray;

	my $total_count = 0;
	for(my $i = 0; $i < @array; $i++){
		my ($target1, $target2) = @{$array[$i]};

		if(abs($source1 - $target1) < 6000 and abs($source2 - $target2) < 6000){
			$total_count++;
			$already_defined{"$target1\t$target2"} = 1;
		}
	}
	return $total_count;
}